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4—* . Abstract. Prior to the recent development of symplectic integrators, 

• ' the time-stepping operator e''' ^^^ was routinely decomposed into a sum 

of products of e*" and e'^^ in the study of hyperbolic partial differential 
equations. In the context of solving Hamiltonian dynamics, we show that 
such a decomposition give rises to both even and odd order Runge-Kutta 
and Nystrom integrators. By use of Suzuki's forward-time derivative op- 
erator to enforce the time-ordered exponential, we show that the same 
f^ \ decomposition can be used to solve non- autonomous equations. In partic- 

fS) ■ ular, odd order algorithms are derived on the basis of a highly non-trivial 

CNj ' time-asymmetric kernel. Such an operator approach provides a general 

and unified basis for understanding structure non- preserving algorithms 
(<— ^ ■ and is especially useful in deriving very high-order algorithms via an- 

("^ ' alytical extrapolations. In this work, algorithms up to the 100th order 

are tested by integrating the ground state wave function of the hydrogen 
atom. For such a singular Coulomb problem, the multi-product expan- 
sion showed uniform convergence and is free of poles usually associated 
^^ ■ with structure-preserving methods. Other examples are also discussed. 

.^, 
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1 Introduction 

In the course of devising numerical algorithms for solving the prototype linear 
hyperbolic equation 



dtu = Aux + Buy, u{0) = uo, (1) 

where A and B arc non-commuting matrices, Strang[34] proposed two second- 
order algorithms corresponding to approximating 

Tih) = c'^('4+B) (2) 

either as 

S(h) = i (c'^^e"^ + c'^^c"^) (3) 

or as 

Following up on Strang's work, Burstein and Mirin[8] suggested that Strang's 
approximations can be generalized to higher orders in the form of a multi-product 
expansion (MPE), 

^HA+B) = ^ cfe [] e'^-'^^e"-''^ (5) 

k i 

and gave two third-order approximations 

nw^l(l^isi^^sM)-lsw (6) 

and 

Bab (h) = ^e^^^/3)A^W3)B^(2h/3)A^(h/3)B _ l^hA^hB^ (7) 
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They credited J. Dunn for finding the decomposition D{h) and noted that the 
weights Cfe are no longer positive beyond second order. Thus the stability of the 
entire algorithm can no longer be inferred from the stability of each component 
product. 

Since (3), (4), (6) and (7) are approximations for the exponential of two 
general operators, they can be applied to problems unrelated to solving hyper- 
bolic partial differential equations. For example, the evolution of any dynamical 
variable w(q, p) (including q and p themselves) is given by the Poisson bracket, 

„ , , /du dH du dH\ ,,„.,. 



For a separable Hamiltonian, 



i^(P,q)-|- + ^(q), (9) 



A and B are Lie operators, or vector fields 



-'-^■^ ^^^"^i; W 



where we have abbreviated v = p/?7i and a(q) = —\7V{q)/m. The exponential 
operators e'''^ and e'^^ are then just shift operators, with S{h) giving the second- 
order Runge-Kutta integrator 



q = qo + /ivo + -/i^a(qo) = qi 



vo+2 



a(qo) + a(qo + /ivq) 
and SAsih), the symplectic Verlet or leap-frog algorithm 

q = qi 

V = vo + - a(qo) + a(q) 
More interestingly, Dunn's decomposition D{h) gives 



q = qo + /ivo + — 
6 



h 
a(qo) -h2a(qo + -vq] 
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a(qo) + 4a(qo + -vq) + 2a(qi) - a(^qi - -h a(qo) 



Since 



2a(qi) - a( qi - -h'^a{qQ)j = a(^qi + -h'^a.{cio) 



0{h^ 



it remains correct to third order to write 



Vo 



h\ 
6 



a(qo) + 4a(qo + ^vq) + a( qo + /ivq + h'^a.{c\f) 



(11) 
(12) 

(13) 
(14) 

(15) 

(16) 

(17) 
(18) 



One recognizes that (15) and (18) as precisely Kutta's third order algorithm[25] 
for solving a second-order differential equation. Burstein and Mirin's approxi- 
mation BAsih) directly gives, without any change, 



q = qo + hvQ 



h^ 



a(qo)+a(q2/3) 



with 



Vo + - a(qo) + 3a(q2/3j 

2 2 

q2/3 = qo + ^hvo + -/i2a(qo), 



(19) 
(20) 

(21) 



which is Nystrom's third order algorithm requiring only two force-evaluations[2, 
31]. Since Burstein and Mirin's approximation is not symmetric, BsAih) pro- 
duces a different algorithm 



q = qo -h ft-vo + yai/g 



Vo 



3ai 



/3 



1/ 
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a(qo + hvo + ^/i^ai/g) - -a(qo + /ivq) 



(22) 

(23) 



where ai/3 == a(qo + /ivo/3). Again, since 

3 4 1 2 

-a(qo + /ivo + -/i2ai/3)--a(qo + /ivo) = a(qo + /ivo + -/t^ai/g) + 0(/i'*), (24) 



(23) can be rewritten as 



V = Vo 



2 1 

3ai/3 +a(qo + /1V0 + -/i^ai/3)J. (25) 



Eqs.(22) and (25) is a new third order algorithm with two force-evaluations but 
without evaluating the force at the starting position. More recently, Ref.[13] has 
shown that Nystrom's four-order algorithm[2] with three force-evaluations and 
Albrecht's six-order algorithni[l] with five-force evaluations can all be derived 
from operator expansions of the form (5). 

Just as symplectic integrators [30, 19] can be derived from a single product 
splitting. 



^h(A+B) ^ 



[]e"''^V'''^, (26) 



these examples clearly show that the multi-product splitting (5) is the funda- 
mental basis for deriving non-symplectic, Nystrom type algorithms. (These are 
not fully Runge-Kutta algorithms, because the operator B in (10) would not 
be a simple shift operator if a(q) becomes dependent on v. On the other hand, 
Nystrom type algorithms are all that are necessary for the study of most Hamil- 
tonian systems.) As illustrated above, one goal of this work is to show that all 
traditional results on Nystrom integrators can be much more simply derived and 
understood on the basis of multi- product splitting. In fact, we have the following 
theorem 

Theorem 1. Every decomposition of c^^ ^^' in the form of 



A^ 



Ck 



YJ^^a.^hA^b.^hB ^ ^h(A+B) ^ C)(/,"+l)^ (27) 



where A and B are non-commmuting operators, with real coefficients {cfe, a^i, hki\ 
and finite indices k and i, produces a nth-order Nystrom integrator. 

(Note that the order n of the integrator is defined with respect to the error in 
approximating the operator {A-\-B) and therefore the error in the time-stepping 
operator is one order higher.) The resulting integrator, however, may not be opti- 
mal. As illustrated above, at low orders, some force evaluations can be combined 
without affecting the order of the integrator. However, such a force consolidation 
is increasely unlikely at higher orders. This theorem produces, both the tradi- 
tional Nystrom integrators where the force is always evaluated initially, and 
non-FASL (First as Last) integrators where the force is never evaluated initially, 
as in (22) and (25). 

The advantage of a single product splitting is that the resulting algorithms 
are structure-preserving, such as being symplectic, unitary, or remain within 



the group manifold. However, single product splittings beyond the second-order 
requires exponentially growing number of operators with unavoidable negative 
coefficients [33, 35] and cannot be applied to time- irreversible or semi-group prob- 
lems. Even for time-reversible systems where negative time steps are not a prob- 
lem, the exponential growth on the number of force evaluations renders high 
order symplectic integrators difficult to derived and expensive to use. For exam- 
ple, it has been found empirically that symplectic algorithms of orders 4, 6, 8 
and 10, required a minimum of 3, 7, 15 and 31 force-evaluations respectively [13]. 
Here, we show that analytically extrapolated algorithms of odd orders 3, 5, 7, 
9 only requires 2, 4, 7, 11 force-evaluations and algorithms of even orders 4, 
6, 8, 10, only require 3, 5, 10, 15 force evaluations. Thus at the tenth order, 
an extrapolated MPE integrator, only requires half the computational effort of 
a symplectic integrator. Or, for 28 force-evaluations, one can use a 14th order 
MPE integrator instead. This is a great advantage in many practical calculations 
where long term accuracy and structure preserving is not an issue. The advan- 
tage is greater still beyond the tenth order, where no symplectic integrators and 
very few RKN algorithms are known. Here, we demonstrated the working of 
MPE algorithms up to the 100th order. 

By use of Suzuki [36] method of implementing the time-ordered exponential, 
this work shows that the multi-product expansion (27) can be easily adopted to 
solve the non-autonomous equation 

dtY{t) = A{t)Y{t), y(0) = Yo. (28) 

In even-order cases, this method reproduces Gragg's[21] classical result in just 
a few lines. In odd-order cases, this method demonstrates a highly non-trivial 
extrapolation of a time- asymmetric kernel, which has never been anticipated 
before. Finally, we show that the multi-product expansion (27) converges uni- 
formly, in contrast to structure-preserving methods, such as the Magnus expan- 
sion, which generally has a finite radius of convergence. The convergence of (27) 
is verified in various analytical and numerical examples, up to the 100th order. 
The paper is outlined as follows: In Section 2, we derive key results of MPE, 
including the extrapolation of odd-order algorithms. In Section 3, we show how 
Suzuki's method can be used to transcribe any splitting scheme for solving non- 
autonomous equations. In Section 4, we present an error and convergence analysis 
of the multi-product splitting based on extrapolation. Numerical examples and 
comparison to the Magnus expansion arc given in Section 5. In Section 6, we 
briefly summarize our results. 

2 Multi-product decomposition 

The multi-product decomposition (5) is obviously more complicated than the sin- 
gle product splitting (26). Fortunately, nineteen years after Burstein and Mirin, 
Sheng[33] proved their observation that beyond second-order, a^i, bki and Ck can- 
not all be positive. This negative result, surprisingly, can be used to completely 
determine Uki, bki and Ck to all orders. This is because for general applications, 



including solving time-irreversible problems, one must have aki and bki positive. 
Therefore every single product in (5) can at most be second-order[33,35]. But 
such a product is easy to construct, because every left-right symmetric single 
product is second-order. Let Ts(h) be such a product with '^^aki = 1 and 
X^i ^ki = Ij then Ts{h) is time-symmetric by construction. 



rs{-h)Ts{h) = 1, 

implying that it has only odd powers of h 

Ts{h) = exp(/i(A + B) + h^Es + h^E^ + ■ 



) 



(29) 



(30) 



and therefore correct to second-order. (The error terms Ei are nested commuta- 
tors oi A and B depending on the specific form of Ts-) This immediately suggests 
that the fcth power of Ts at step size h/k must have the form 

Ti{h/k) = exp(/i(A + B) + k-^h^Ea + k-^h^E^ 



+ •••), 



(31) 



and can serve as a basis for the multi-production expansion (5). The simplest 
such symmetric product is 



r2{h)^SAB{h) or T2{h) = SsAih). 
If one naively assumes that 

T2{h) = e''(^+^) +Ch'' + Dh'' + --- , 
then a Richardson extrapolation would only give 



1 



^^^-^[k%\h/k)-T2{h)]=C 



_ MA+B) 



0{h^), 



(32) 
(33) 

(34) 



a third-order [32] algorithm. However, because the error structure oi T2(h/k) is 
actually given by (31), one has 



r2'=(Vfc)-e''(^+'''- 



k-^h'^Ea + lk-^h'^[{A + B)E3 + E3{A + B)] + 0{h'^), (35) 



and both the third and fourth order errors can be eliminated simultaneously, 
yielding a fourth-order algorithm. Similarly, the leading 2n + l and 2n + 2 order 
errors are multiplied by /c^^" and can be eliminated at the same time. Thus for 
a given set of n whole numbers {ki} one can have a 2nth-order approximation 



^HA+B) 



Yl ^^'^2' 



i=l 



0{h 



2n+l-\ 



(36) 



provided that Ci satisfy the simple Vandermonde equation: 



/ 1 



1 



1 



,-2("-l) L-2("-l) 



1 

/r-2 



il,-2("-l) 



\ 


C2 
C3 


= 


/ 1 \ 





) 


UJ 




Uy 



(37) 



Surprisely, this equation has closed form solutions [13] for all n 

The natural sequence {ki} = {1,2,3 ...n} produces a 2nth-order algorithm with 
the minimum n{n + l)/2 evaluations of T2{h). For orders four to ten, one has 
explicitly: 



%{h) 



-ir,WH-i7?(i) (39) 

x,,.,.lr,.,-|r/(i).|7,»(i) (40, 

16384 , /fe\ 390625 , fh\ 

As shown in Ref.[13], 71 (/i) reproduces Nystrom's fourth-order algorithm with 
three force-evaluations and Te{h) yielded a new sixth-order Nystrom type algo- 
rithm with five force-evaluations. 

Remark 1. It is easy to show that the Verlet algorithm (13) and (14) correspond- 
ing to SAsih) produces the same trajectory as Stomer's second order scheme 

qi - 2qo + q_i = h^a{qo). (43) 

However, it is extremely difficult to deduce from (43) that the underlying error 
structure is basically (30) and allows for a /i^-extrapolation. This is the great 
achievment of Gragg[21]. Nevertheless, the power of the present operator ap- 
proach is that we can reproduce his results in a few lines. The error structure 
here, (30), is a simple consequence of the symmetric character of the product, 
and allows us to bypass Gragg's lengthy proof on the asymptotic errors of (43). 
Moreover, this ^.^-extrapolation can be applied to any Ts{h), not necessarily re- 
stricted to (43). For example, the use of SsAih) produces an entirely different 
sequence of extrapolations [13], distinct from from that based on (43). 

Remark 2. In the original work of Gragg, the use of (43) as the basis for his 
extrapolation is a matter of default; it is a well-known second-order solution. 
Here, in extrapolating operators, the use of SAsih) or SsAih) is for the specific 
purpose that they can be applied to time-irreversible problems. While all pos- 
itive time steps algorithms are possible in the fourth-order [3 7, 10] by including 
the operator [B, [A, B]], MPE is currently the only way of producing sixth and 



higher-order algorithms in solving the imaginary time Schrodinger equation[14] 
and in doing Path- Integral Monte Carlo simulations [41]. The fact that MPE is no 
longer norm preserving nor even strictly positive, does not affect the higher order 
convergences in these applications. These non-structure preserving elements are 
within the error noise of the algorithm. MPE is less useful in solving the real 
time Schroding equation where, unitarity is of critical importance. 

Remark 3. The explicit coefficient Ci coincide with the diagonal elements of the 
Richardson- Aitken-Neville extrapolation[24] table. This is not surprising, since 
they are coefficients of extrapolation. As shown in [13], q — Li{0), where Li{x) 
are the Lagrange interpolating polynomials with interpolation points Xi = k~ . 
What is novel here is that Ci is known analytically and a simple routine calling 
it repeatedly to execute 72 (/i) will generate an arbitrary even order algorithm 
without any table construction. The resulting algorithm is extremely portable 
and compact and can serve as a benchmark by which all integrators of the same 
order can be compared. In Ref. [13], the only algorithm that have outperformed 
MPE is Dormand and Prince's [15] 12th-order integrator as given in Ref. [7]. 

Having the explicit solutions Ci now suggests new ways of solving old prob- 
lems. Suppose one wishes to integrate the system to time t. One may begin by 
using a second-order algorthm and iterate it m time at time step h = t/m, 

T2,m{h)^Tr{tlm). (44) 

Every position on the trajectory will then be correct to second order in h. How- 
ever, if one were only interested in the final position at time i, then one can 
correct this final position to fourth order by simply computing one more T2 (t) 
and modify (44) via 

%Mh) = -^^mt/m) -^^T2(t), (45) 

?7l^ — 1^ 771^ — 1^ 

or correct it to sixth-order via 

m*Tr{t/m) 247?(i/2) 1^X2(0 



Tlo,m{h) 



(to2 _ 12)(^2 _ 22) (22-l2)(22-m2) (P - 22) (l2 - ^2) ' 

(46) 
and so on, to any even order. The expansion coefficients are given by {ki\ equal 
to {to, 1}, {to, 2,1}, {to, 3, 2,1} etc.. This is similar to the idea of process al- 
gorithms[3], but much, much simpler. The processor for correcting (44) beyond 
the fourth-order can be quite complex if the entire algorithm were to remain 
symplectic. Here, for Nystrom integrators, the extrapolation coefficient is known 
to all even orders. Alternatively, one can view the above as correcting every TO.th 
step of the basic algorithm T2{t/m) over a short time interval oft. Thus knowing 
Ci allows great flexibility is designing algorithms that run the gamut from being 
correct to arbitrary high order at every time-step, every other time-step, every 
third time-step, etc., to only at the final time step. With MPE, one can easily 
produce versatile adaptive algorithms by varying both the time step size h and 
the order of the algorithm. 



Remark 4- Since MPE is an extrapolation, it is expected to be more proned to 
round-ofF errors. Thus if n is too large in (45), the second term maybe too small 
and the correction is lost to round-off errors. However, as seen in (17) and (24), 
the required substractions are sometime well-defined and the the round-off errors 
are within the error noise of the algorithm. As will be shown in Section 5, the 
round-off errors are sometime less severe than expected. 

Remark 5. The idea of extrapolating symplectic algorithms has been considered 
by previously by Blanes, Casas and Ros[6] and Chan and Murua[9]. They studied 
the case of extrapolating an 2n-order symplectic integrator. They did not obtain 
analytical forms for their expansion coefficients but noted that extrapolating 
a 2n-order symplectic integrator will preserve the symplectic character of the 
algorithm to order 4ri + 1. While this is more general, such an extrapolation 
cannot be applied to time- irreversible systems for n > 1. 

Finally, we note that 

i—l ^ ' 

In principle, for any countable sets of {fci}, we have achieved an exact decompo- 
sition, with known coefficients. This is in contrast to the structure-preserving, 
but impractical Zassenhaus formula. 

The above derivation of even-order algorithms, is at most an elaboration 
on Gragg's seminal work. Below, we will derive arbitrary odd-order Nystrom 
algorithms which have not been anticipated in any classical study. Since 

Tiih) = e'^'^e''-^ = exp[/i(A + B) + h^F2 + h^Fs + h'^F^ + •••], (48) 

contain errors of all orders {{Fi\ are nested commutators of the usual Baker- 
Campbell-Hausdorff formula), extrapolations based on Ti{h/k) will not yield 
a /i^-order scheme. However, there is a /i^-order basis hidden in Burstein and 
Mirin's original decomposition (7). The following basis for n = 1, 2, 3 . . . 

U,^(h) = e5^^(e^^^e^^-4)"-ie5^^ (49) 

has the remarkable property that it effectively behaves as if 

Un{h) == cxp[h(A + B)+ x-\h^F2 + h^Fa) + x'^h'^Fi + h'^F^) + • • • ] (50) 

where x — {2n — 1). (By effectively we mean that Un{h) actually has the form 

Un(h) ^ exp[h{A + B)+ x-^{h^F2 + h^F^) 

+ (a;-2 _ x-^)h''F'i + x-^h^Fi + h^'F^) + • • • ] (51) 

where F!^ are additional commutators not present in (48). However, this is es- 
sentially (50) with altered Fi but without changing the crucial power pattern of 
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x^'^''.) In this case, (50) (as well as (51)) can be extrapolated similarly as in the 
even order case, 

n 

^HA+B) ^ J2 £^U,{h) + 0{h^^), (52) 

where Ci satisfies the same Vandermonde equation (37), with the same solution 
(38), but with {ki} consists of only odd whole numbers. The first few odd or- 
der decompositions corresponding to {ki} being {1,3}, {1,3,5}, {1,3,5,7} and 
{1,3,5,7,9} are: 

Ts{h)^^^U,{h) + ^U2{h) (53) 

1 R^ ft'?^ 

%{h) = —U,{h) —U,{h) + —^U,{h) (54) 

1 72Q 1 5625 1 1 7fi4Q 

^^ ^ 9216 ^ ^5120 ^ ' 9216 ^^ ' 46080 ^ ' ^ ' 

^r,. 1 ., ,,. 729 ,, ,,, 390625,, ,, , 

1474560 ^ ' 1146880 ^^ ' ^ ' 

The splitting Taih) explains the original form of Burstcin and Mirin's de- 
composition and Nystrom's third-order algorithm. The splitting T^ih) again 
produces, without any tinkering, Nystrom's fifth-order integrators [2] with four 
force-evaluations : 



1 + 
h 



q = qo + /ivo + — |^23ao + 75a2/5 - 27a2/3 + 25a4/5j (57) 



V = Vn H 

192 



23ao + 125a2/5 - 81a2/3 + 125a4/5j , (58) 

where we have denoted aj/j, = a(qj/fe) with 

2 2 

q2/5 = qo + -r/lVo + ^h^^o 

4 4 o 

q4/5 = qo + -r/ivo + — /i (ao + a2/5) (59) 

and where q2/3 has been given earlier in (21). (Interchange of A <-)• _B in T5{h) 
will also yield a fifth-order algorithm, but since the final force-evaluations can 
only be combined as in (24) to order 0{h*), such a force consolidation cannot 
be used for a fifth-order algorithm. The algorithm will then require six force- 
evaluations, which is undesirable. We shall therefore ignore this alternative case 
from now on.) With three more force-evaluations at 



q2/7 = 


2, 2 ^2 
= qo + i=hvo + —h ao 
7 49 


q4/7 = 


= qo + y/ivo + ^^^(ao + a2/7) 


qe/T = 


6 2 
= qo + y/ivo + ^^^(3ao + 4a2/7 + 2a4/7 



(60) 
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Trih) produces the following seventh-order algorithm with seven force-evaluations, 
which has never been derived before, 

q = qo + hvo 



23040 



h 

vo 



23040 



1682ao + 729a2/3 - 3125(3a2/5 + a4/5) 

+2401(5a2/7 + 3a4/7 + a6/7)] (61) 
1682ao + 2167a2/3 - 15625(a2/5 + 3.^/,^) 

+16807(a2/7 + a4/7 + a6/7)l • (62) 



These analytical derivations are of course unnecessary in practical applications. 
As in the even-order case, both the coefficients Ck and the algorithm correspond- 
ing to hin{h) can be called repeatedly to generate any odd-order integrators. 

Since each lAn{h) requires n force evaluation, but have the initial force in 
common, each (2n — 1) order algorithm requires ^n{n — 1) -|- 1 force-evaluations. 
Thus for odd-orders 3, 5, 7, 9, the number of force-evaluation required are 2, 4, 7, 
11. As alluded to earlier, for even-order 4, 6, 8, 10, the number of force-evaluation 
required are 3, 5, 10, 15. These sequences of extrapolated algorithms therefore 
provide a natural explanation for the order barrier in Nystrom algorithms. For 
order p < 7, the number of force-evaluation can he p — 1, but for p > 7, the 
number of force-evaluation must be greater than p. 

Remark 6. In general we have the following order notation for the even and odd 
algorithms: 

— The order of the even algorithm is 2n, its decomposition error is 2n + 1. 

— The order of the odd algorithm is 2n — 1, its decomposition error is 2n. 

3 Solving non-autonomous equations 

The solution to the non- autonomous equation (28) can be formally written as 

Y{t + h)=r(cxp I A{s)ds)Y{t), (63) 

aside from the conventional expansion 

r-t+h r-t+h r-t+h r-Si 

rf exp / A{s)ds) ^1+ A{si)dsi+ dsi ds2A{si)A{s2) + - ■ ■ , 

(64) 
the time-ordered exponential can also be interpreted more intuitively as 

r(exp / A{s)ds) = Imi rfe^2:r=iA(t+«i-)\ (65) 

= lim et^(*+'')...e^^(*+^)e*-4(*+^). (66) 
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The time-ordering is trivially accomplished in going from (65) to (66). To enforce 
latter, Suzuki[36] introduces the forward time derivative operator, also called 
super-operator: 

such that for any two time-dependent functions F{t) and G{t), 

F{t)e''^G{t) = F{t + h)G{t). (68) 

If F(t) == 1, we have 

le''^G(t) = c^'^Git) ^ G{t). (69) 

Trotter's formula then gives 

eMHMt) + D)] = lim fe^^We-^r, 

- lim e^^(*+")...e^-4(*+^)e^^(*+^), (70) 

n— 7-00 

where property (69) has been applied repeatedly and accumulatively. Comparing 
(66) with (70) yields Suzuki's decomposition of the time-ordered exponential[36] 

rfexp /" A{s)ds) ^ cxp[h{A{t) + D)]. (71) 

Thus time-ordering can be achieve by splitting an additional operator D. This is 
extremely useful and transforms any existing splitting algorithms into integrators 
of non-autonomous equations. For example, one has the following symmetric 
splitting 

T2{h) = c^'^D^hA(t)^^hD ^ ^hAit+^h)^ (72) 

which is the second-order mid-point approximation. Every occurrence of the 
operator e'^^^^ ^ from right to left, updates the current time t to t + dih. If t is 
the time at the start of the algorithm, then after the first occurrence of e^^^, 
time is t + i/i. After the second e^''^, time is t -t- ft-. Thus the leftmost e^''^ is 
not without effect, it correctly updates the time for the next iteration. Thus the 
iterations of 72 (ft) implicitly imply 



(73) 



For the odd-order basis, we have 

Z^i(ft)=e'*^e'*'4(t)^g/>A(t) 

Usih) = ei:hA{t+l:h)^%hA{t+lh)^lhAit) 



(74) 
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Remark 7. The recent work by Wiebe et al.[38] suggests that Suzuki's decompo- 
sition (71) only holds li A{t) is sufSciently smooth. In cases where the derivatives 
of A{t) cease to exist, high-order integrators based on (71) maybe degraded to 
lower orders. 

For A{t) = T + V{t), since [D,T] ~ 0, the second-order algorithm can be 
obtained as 

T2{h) = c^HT+D)^hV(t)^^h{T+D) 

= e^''^e''^(*+''/^)e^''^. (75) 

For odd order algorithms, we now have the following sequence of basis product 

Z^i(/j) = e'^^e'^^(*) 

U2{h) = e\hT^lhV(t+lh)^lhT^^hV{t) 

Ui{h) = eJ'*^eJ''^(*+^'')e*''^e^''^(*+j'')eJ'''^ei''^(*) 

(76) 

While any power of 72 (^) is timc-symmertic, each Un{h) is time asymmetric, 

Un{-h)Un{h) ^ 1. (77) 

4 Errors and convergence of the Multi-product expansion 

While extrapolation methods are well-known in the study of differential equa- 
tions, there is a virtually no work done in the context of operators. Here, we ex- 
tend the method of extrapolation to the decomposition of two operators, which 
is the basis of the MPE method. Working at the operator, rather than at the 
solution level, allows the extrapolation method be widely applied to many time- 
dependent equations. In particular, we will use the constructive details in[13] to 
prove convergence results for the multi-product expansion. While this work is 
restricted to exponential splitting, our proof of convergence based on the general 
framework of [22]. 

4.1 Analysis of the even-order kernel 72 

We will assume that at sufficient small h, the Strang splitting is bounded as 
follow: 

\\r2{h)\\ = \\exp{^hD)cxp{hA{t))exp{-hD)\\ < cxp{cLjh), (78) 

with c only depend on the coefficients of the method, see the work of convergence 
analysis on this splitting by Janke and Lubich [27]. We can then derive the 
following convergence results for the multi-product expansion. 
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Theorem 2. For the numerical solution of (28), we consider the MPE algo- 
rithm (36) of order 2n. Further we assume the error estimate in equation (78), 
then we have the following convergence result: 

II (5" - eMmh{A{t) + D))) uo\\ < CO{h^^+'),mh < te„d, (79) 

where S = "^"^i CiT2 ^ {-r-) and C is to be chosen uniformly on bounded time 
intervals and independent of m and h for sufficient small h. 

Proof. We apply the telescopic identity and obtain: 

(S"" - exp{mh{A{t) + D))) uq = (80) 

J2 S"'-''-\S - cxp{h{A{t) + D))) cxp{i^h{A{t) + D))uo. (81) 

where S = J:tlC^Ti'iI-) 

We apply the error estimate in (78) to obtain the stability requirement: 

\\Y.c,r^^{-)\\<exj>{cuh). (82) 

Assuming the consistency of 

" h 

llE^^7-2''(r)-'^^PW^ + ^))ll ^C'0(/i2"+i) (83) 

we have the following error bound: 

II (5" - exp{mh{A{t) + D))) uo\\ < CO{h^''+^),mh < t^nd, (84) 

The consistency of the error bound is derived in the following theorem. 

Theorem 3. For the numerical solution of (28), we have the following consis- 
tency: 

\\J2^^'^2'iTr)-^MhiA + D))\\<CO{h'^^+'). (85) 

Proof. Based on the derivation of the coefficients via the Vandermonde equation 
the product is bounded and we have: 

E ^kT.H^) = E ^* (^^P((^ + ^)'^) - (k-^h^Es + k-^h^E, + ...)) ,(86) 

n / n \ 

= ^ c J exp((yl + D)h) - ^ k-^'h^'+'E2^+i , 

(n n \ 

exp((A + D)h) - E C'^ E k-^'h^'^'E2^+l , 
fe=l i=l / 

= 0(/i2"+i), (87) 
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where the coefficients are given in (38). 

Lemma 1. We assume \\A{t)\\ to be bounded in the interval t G (0,ie„(i)- Then 
72 is non- singular for sufficient small h. 

Proof. We use our assumption ||A(t)|| is to be bounded in the interval < i < 

So we can find ||A(i)|| < C for < i < ignd, where C G IR+ a bound of 
operator A{t) independent of t. 

Therefore 72 is always non-singular for sufficiently small h. 

Remark 8. Based on these results the kernel 72 is also uniform convergent. 

The same argument can be used by applying to MPE formula, while all 
operators are convergent, the sum of all is also bounded and convergent, see [16] 
and [18]. 

Remark 9. For higher kernels, e.g. 4th order, there exists also error bounds so 
that uniform convergent results can be derived, see e.g. [20]. Such kernels can 
also be used to the MPE method to acchieve higher order accuracy with uniform 
convergent series. But as we noted earlier, these cannot be applied to time- 
irreversible problems. 

4.2 Analysis of the odd-order kernel Un 

Lemma 2. We will assume that for sufficiently small h, the Burstein and Mirin's 
decomposition is bounded as follow: 

\\Un{h)\\ = lle^^'4(t)(g,|^Dg5|^A(t)y.-i^5^D|| < cxp(cco/^), Vt > 0,(88) 
with c only dependent on the coefficients of the method. 

The proof follows by rewriting equation (88) as a product of the Strang and 
the A-B splitting schemes: 

Proof. Equation (88) can be rewritten as: 

e2;fcr^(*)(e2arr^e^^^(*))"-ie5^^ (89) 

The error bound and underlying convergence analysis for both the Strang and 
the A-B splitting have been previously studied by Janke and Lubich [27]. 

We assume the following derivation of the higher order MPE: 

Assumption 1 We assume the following higher order decomposition, 



^h(A+D) _ 



= ^5,Z^,(/i)+0(/i'"). (90) 



where Ci are derived based on the Vandermonde equation (37) with {ki\ being a 
set of odd whole numbers. 
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We can then derive the following convergence results for the multi-product 
expansion. 

Theorem 4. For the numerical solution of (28), we consider the Assumption 1 
of order 2n — 1 and we apply Lemma 2, then we have a convergence result given 
as: 

\\{S"'-exp(mh{A{t)+D)))uo\\ < CO{h^''),mh < t^nd, (91) 

with n = 1, 2, 3, . . ., and where S = X]"=i CiUi{h) and C is to he chosen uniformly 
on bounded time intervals and independent of m and h for sufficient small h. 

Proof. The same proof ideas can be followed after the proof of Theorem 2. 
The consistency of the error bound is derived in the following theorem. 

Theorem 5. For the numerical solution of (28), we have the following consis- 
tency: 

n 

II ^ Q U,{h) - exp{h{A + D))\\< CO{h^n- (92) 

i=l 

Proof. The same proof ideas can be followed after the proof of Theorem 3. 

Remark 10. The same proof idea can be used to generalise the higher order 
schemes. 



5 Analytical and numerical verifications 

In this section, we seek to verify and assess the convergence of both the even 
and odd order MPE algorithms. For a single product splitting, there are no 
known splittings that are exact in the limit of large number of operators. Even 
in the case of the Zassenhaus formula, it is non-trivial to compute the higher 
order products, not to mention evaluating them. For this purpose, wc turn to 
the much studied Magnus expansion, where the exact limit can be computed in 
simple cases. 

The Magnus expansion[5] solves (28) in the form 

oo 

Y{t) = exp(i7(t))r(0), f2{t) ^ ^ I2„(t) (93) 

where the first few terms are 

Qi{t) = I dtiAi 
Jo 

^2(t)-i / dti [ dt2[Ai,A2] 

^ Jo Jo 



n=l 



^3(t) = ^ / dti r dt2 r dt^{[Ai,[A2,A^] + [[Ai,A2],A^]) 

D Jo Jo Jo 



(94) 
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with An = A{tn). In practice, it is more useful to define the nth order Magnus 
operator 

I2["l(t) = i7(i)+0(i"+i) (95) 

such that 

Y{t) = exp[r2["l(t)]r(0) + 0(r+i). (96) 

Thus the second-order Magnus operator is 







= tA[^t]+ 0{f) (97) 



and a fourth-order Magnus operator [5] is 



l2W(t) ^ it(Ai + A2) -c3i2[Ai,A2] (98) 



C3 - ^. (99) 



where Ai = A(cii) 


,^2 


= A{C2 


t) and 






Cl 


1 

^ 2 


V3 
6 ' 


C2 


1 

= 2 + 


V3 
6 


For the ubiquitous 


case 


of 














A(i) = 


:T + F(i), 


one has 














e 


n^^\t) ^ 


= gt[T+y(t/2)] 
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(100) 



^gitTgty(t/2)gi*T^O(t3) (101) 

and 

ef2W(t) ^ gC3t(y2-Vi)gt(T+i(yi+y2))g-c3t(y2-Vi) _j_ ^,^^5) ^^q2) 

where 

Vi = F(cii), V2 - y(c2i). (103) 

The Magnus expansion (96) is automatically sturcture-preserving because it is 
a single exponential operator approximation. However, since one must further 
split i?!"! into computable parts, the expansion is as complex, if not more so, 
than a single product splitting. In the following, the comparison is not strictly 
equitable, because the MPE is not structure-preserving. Neveetheless it is useful 
to know that, perhaps for that reason, MPE can be uniformly convergent. 

5.1 The non-singular matrix case 

To assess the convergence of the Multi-product expansion with that of the Mag- 
nus series, consider the well known example[29] of 

A{t)=(l^\. (104) 
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The exact solution to (28) with Y{0) = / is 

■32* /(t) 



Yit) = 







with 



/(i) = _e-*(e^*_l-3i) 



2 



f* 



i-^ 



i" 



t' 



31t« 



t« 



13^1 



13^1 



8 60 80 420 40320 6720 403200 178200 
For the Magnus expansion, one has the series 



(105) 



(106) 
{107) 



f2{t) 



2i g{t) 
-t 



with, up to the 10th order, 

git) 



9 



2 4 80 1120 



t' + 



81 



44800 



e + ' 



l-3t) 
3(e3* - 1) ■ 

Exponentiating (108) yields (105) with 

1 1 1 



1/^31 



fit) = te-'(e' 



^ * 6 12 80 1120 



3 . 27 



44800 



t' + 



^te-'(c^*-l) -- 



1 



1 



9i 3(e3* - 1) 



(108) 

(109) 
(110) 

(111) 
(112) 



Whereas the exact solution (106) is an entire function of i, the Magnus series 
(109) and (111) only converge for \t\ < |7r due to the pole at t = |7rz. The 
Magnus series (111) is plot in Fig.l as blue lines. The pole at |i| = |7r ss 2 is 
clearly visible. 

For the even order multi-product expansion, from (72), by setting h = t and 
i = 0, we have 



T2{t) =exp t(^ 2 
and we compute T2{t) according to (73) as 
7;2(t/2)=exp 
with 



,2t 



t (2 \t 



2 VO -1 
1 



./2(t) 
0-* 



C^2 / 



/2(t)-^te-*(e^*-l) 
6 



(113) 

(114) 
(115) 
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20 



15 - 



10 - 



5 - 



- 




Fig. 1. The black line is the exact result (106). The dotted blue lines are the Magnus 
fourth to tenth order results (111), which diverge from the exact result beyond i > 2. 
The solid red lines are the multi-product expansions. The dashed-purple line is their 
common second order result. 



This is identical to first term of the Magnus series (111) and is an entire func- 
tion of t. Since higher order MPE uses only powers of 72, higher order MPE 
approximations are also entire functions of t. Thus up to the 10th order, one 
finds 

.Uit)=t^-'l-^ + ^] (116) 



18 



9 



hit) = te"* 1 e^* -He* e^*/2 

■'^^ ' ^ 360 40^ ^ 45 



(117) 



fsit)=te- 



151e3* - 2369 256 



81 



7560 



+ ^(e«*/* + e^*/*) - ---(e^' + e*) + 



945 



280 



104 
315^ 



,3t/2 



(118) 



^ ,, _t/15619e3*- 347261 78125 , i2t/5 9t/5 6t/5 3t/^i^ 
•'^°^ ^ \ 1088640 217728^ -r -r -r ; 



4096, 



729 



4192 



:(e9*/4 + e3*/4) + ^±^(e2*+e*) 
8505^ ' 4480^ ^ 8505 



,3t/2 



(119) 
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These even order approximations are plotted as red lines in Fig.l. The conver- 
gence is uniform for all t. 

When expanded, the above yields 

" t3 



t 



h{t) = ^ + - + 



hit) - I 



hit) = IT + ^ + 7^ + 



hit) - y 



/io(t) 



2 



5t^ 
192 

60 

60 

60 



80 
80 
80 



384 
420 
420 



31t" 

40320 

31t8 



1307f^ 

8601600 + '" 
t9 13^10 



13099^1 



7^120) 



40320 6720 403200 232243200' 
and agree with the exact solution to the claimed order. Similarly, the m-step 
extrapolated algorithms 72, m, 74,™ , etc., are also correct up to the claimed order. 
For odd orders, by again setting h — t and t — 0, the basis defined in (74) 
now reads 



Uiit) 

U2{t) 



cxp 



cxp 



e2* 
e-* 



^ 3' 

-1 

W - e" 



exp 



1 



-t 



2 
-1 



(121) 



and the MPE (53) to (56) give 



hit) 



hit) = y 

t^ 
hit) - ^ 

hit) - % 



12 



t" 



t^ 

60 


+ 


llt^ 
1000 


+ •• 












t^ 
60 


+ 


t^ 

80 + 


420 


+ 


18299i^ 
24696000 








t^ 


+ 


t6 

80 + 


420 


+ 


31f« 
40320 + 


t^ 


+ 


1577*10 


60 


3720 


49392000 



.(122) 



Results (120) and (122) constitute an analytical verification of the even and odd 
order MPE (39)-(42) and (53)-(56). 

5.2 The singular matrix case 

Consider the radial Schrodinger equation 



Qj.2 



!{r,E)u{r) 



(123) 
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where 

f{r,E)^2V{r)-2E+^-^i^. (124) 

By relabeling r ^ t and u{r) — > q{t), (123) can be viewed as harmonic oscillator 
with a time dependent spring constant 

kit,E)^-fit,E) (125) 

and Hamiltonian 

H = ^p^ + h{t,E)q\ (126) 

Thus any eigenlunction ol (123) is an exact time-dependent solution of (126). 
For example, the ground state of the hydrogen atom with I ^ 0, E = —1/2 and 

Vir)^-- (127) 

r 

yields the exact solution 

q{t) = texp(— i) 

2 t^ t^ t^ t^ f t« 

= t-r + - 



I 6 24 120 720 5040 
= t-t^ + 0.1667f'' + 0.0417t^ - 0.0083t'^ • • • (128) 

with initial values q{Q) — and p{Q) — 1. Denoting 

the time-dependent harmonic oscillator (126) now corresponds to 



with a singular matrix element 



t 
The second-order midpoint algorithm is 

r2(/i,t)-e5'*^e''^(*+V2)e^"^ 



fit) = (1 - ?). (131) 



hf{t+\h) i + i^fit+m, 



(132) 
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and for q(0) ~ and p{0) ~ 1, (setting t = and h = t), correctly gives the 
second order result, 

q,it)^t+h^f{^t) = t-e + ^t\ (133) 

The even order multi-product expansions (39)-(42) then yield 

qi{t) =t-t^ + 0.3889i^ - O.lllli'' + 0.0104^"^ 
qe{t) =t-t^ + 0.4689^3 - 0.1378i'' + 0.0283*'"^ - 0.0043i^ 
qs{t) =t-t^ + 0.4873^3 - 0.1542i'' + 0.0356t^ - 0.0062^^ • • • 
gio(i) ^1-1"^ + 0.4936i^ - 0.1603i'' + 0.0385t^ - 0.0073i^ • • • (134) 

where we have converted fractions to decimal forms for easier comparison with 
the exact solution (128). One sees that MPE no longer matches the Taylor ex- 
pansion beyond second-order. This is due to the singular nature of the Coulomb 
potential, which makes the problem a challenge to solve. (If one naively makes 
a Taylor expansion about t — starting with q{0) = and p(Q) = 1, then every 
term beyond the initial values would either be divergent or undcfine.) 

Since A(t) is now singular at i = 0, the previous proof of uniform convergence 
no longer holds. Nevertheless, from the exact solution (128), one sees that force 
(or acceleration) 

liinf{t)q{t) = -2 (135) 

remains finite. It seems that this is sufficient for uniform convergence as the 
coefficients of t" do approach their correct value with increasing order. 

For odd order MPE, while each term e^'*/^-'^^*-' of the basis product in (76) 
is singular at t = 0, but because of (135), 

L"r'"""'"(:S)-(i-2v.)' <'=^' 

Interpreting the action of the first operator this way, the basis products of (76) 
then yield, according to the MPE (53)- (56), 

qsit) ^t-f + — -o.nnt* 

/3 

qsit) =t-f + — - 0.1458f'' + 0.0333t^ - 0.0033i^ 

qrit) ^t-f + — - 0.1628f'' + 0.0382t^ - 0.0067^^ • • • 

qgit) =t-t^ + — - 0.1655f'' + 0.0406t^ - 0.0078^^ • • • (137) 

Now 93 (i) is correct to third order, but higher order algorithms are still down- 
graded and only approaches the exact solution asymptotically but uniformly. 

To see this uniform convergence, we show in Fig. 2, how higher order MPE, 
both even and odd, up to the 100th order, compares with the exact solution. The 
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Fig. 2. The uniform convergence of the muhi-product expansion in solving for the 
hydrogen ground state wave function. The black line is the exact ground state wave 
function. The numbers denote the order of the multi-product expansion. The dotted 
blue lines denote results of various fourth-order algorithms. 



calculation is done numerically rather than by evaluating the analytical expres- 
sions such as (134) or (137). The order of the MPE algorithms are indicated by 
numbers. For odd order algorithms, we do not even bother to incorporate (136), 
but just avoid the singularity by starting the algorithm at i = 10~^. Also shovifn 
are some well know fourth-order symplectic algorithm FR (Forest- Ruth[19], 3 
force-evaluations), M (McLachlan[28], 4 force-evaluations), BM (Blanes-Moan[4], 
6 force-evaluations), Mag4 (Magnus integrator, 4 force-evaluations) and 4B[12] 
(a forward symplectic algorithm with sa 2 evaluations). These symplectic in- 
tegrators steadily improves from FR, to M, to Mag4, to BM to 4B. Forward 
algorithm 4B is noteworthy in that it is the only fourth-order algorithm that 
can go around the wave function maximum at t — 1, yielding 



^3 

q4B(i) ^t-t^ + — - 0.1635t^ 



0.0397t^ - O.OOTOt^ + 0.0009t^ • 



(138) 



with the correct third-order coefficient and comparable higher order coefficients 
as the exact solution (128). By contrast, the FR algorithm, which is well know 
to have rather large errors, has the expansion, 



qpRit) ^t~t^ - 0.1942*3 -(- 3.528*'' - 2.415*^ -I- 0.5742*^ - 0.0437i^ 



(139) 
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with terms of the wrong signs beyond t^ . The failure of these fourth-order algo- 
rithms to converge correctly due to the singular nature of the Coulomb potential 
is consistent with the findings of Wiebe et aL[38]. However, their finding does not 
explain why the second-order algorithm can converge correctly and only higher 
order algorithms fail. A deeper understanding of Suzuki's method is necessary 
to resolve this issue. 




Fig. 3. The multi-product expansion results for the hydrogen ground state wave func- 
tion using only double-precision arithematics. The effect of limited precision begins to 
show at order 60. The results for order 65 and 70 are indicated as black and dotted 
blue lines respectively. 



In Fig. 2, results for orders 60, 80 and 100 are computed using quadruple 
precision. If one uses only double-precision, the effect of round-off errors on 
limited precision is as shown in Fig. 3. For this calculation, the round-off errors 
are not very noticeable even at orders as high as 40 or 49, which is rather 
surprising. The round-off errors are noticeable only at orders greater than « 60. 

For non-singular potentials such as the radial harmonic oscillator with 



f{t) ^f~3, 



(140) 
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and exact ground state solution 

+ 2 /r) t t t t t , ^ 

'^W-^^ ' =*-Y + ¥-48 + 384-3840 + ---' ^'''^ 

the multi-product expansion has no problem in reproducing the exact solution 
to the claimed order: 



q6(t) = t-^ 



t^ f5 i^f 



2 8 576 

t^ f 1082t9 



97 (t) = 


= t- 


-J + J- 


48 


98 (t) = 


= t- 


<3 i5 


48 


99 (t) = 


= t- 


t3 i5 


" 48 



385875 
20803*9 
7741440 ^ " ' 
t^ 341f" 



384 1224720 
t^ f5 ^7 ^9 50977i" 

gin(t) =f \ h ^ . (142) 

^ '^ ' 2 8 48 384 193536000 ^ ' 

In this case, the odd order algorithms have the advantage of being correct one 
order higher. 



6 Concluding summary and discussions 

In this work, we have shown that the most general framework for deriving 
Nystrom type algorithms for solving autonomous and non-autonomous equa- 
tions is multi-product splitting. By expanding on a suitable basis of operators, 
the resulting multi-product expansion not only can reproduce conventional ex- 
trapolated integrators of even-order but can also yield new odd-order algorithms. 
By use of Suzuki's rule of incorporating the time-ordered exponential, any multi- 
product splitting algorithm can be adopted for solving explicitly time-dependent 
problems. The analytically know expansion coefficients Ci allow great flexibility in 
designing adaptive algorithms. Unlike structure-preserving methods, such as the 
Magnus expansion, which has a finite radius of convergence, our multi-product 
expansion converges uniformly. Moreover, MPE requires far less operators at 
higher orders than either the Magnus expansion or conventional single-product 
splittings. The general order-condition for multi-product splitting is not known 
and should be developed. In the future we will focus on applying MPE methods 
for solving nonlinear differential equations and time-irreversible or semi-group 
problems. 
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